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dL. INTRODUCTION 


A. BACKGROUND 

Autonomous vehicles suitable for use in modern 
applications require high maneuverability, quick response, and 
robust performance characteristics [Ref. 7]. in -erder to 
maintain accurate track keeping characteristics, a suitable 
combination of path planning, navigation, guidance, and 
control design is needed [Ref. 10]. Sufficient information is 
obtained from charted obstacles and the environment for smooth 
paths to be generated for the vehicle to follow [Ref. 8]. 
Although it is possible to design a nonlinear controller which 
incorporates and utilizes information on the nonlinear dynamic 
properties of the vehicle, as well as the geometric 
nonlinearities of the desired track, the resulting scheme 
tends to be very complex and time consuming [Ref. 3]. The 
alternative is to use separated navigation, guidance, and 
autopilot functions. A certain level of feedback is provided 
at the navigation level through the use of sonars in order to 
replan a path due to uncharted obstacles or changes in mission 
objectives. This operation is event-driven and occurs 
asynchronously, and in many cases it does not affect stability 
and performance of the overall navigation, guidance, and 


control scheme. Based on the provided navigational 


information, the guidance law generates heading commands, 
which are executed by the control law by appropriate use of 
vehicle effectors, such as control surfaces and cross body 
thrusters. However, the time required to process sonar data 
and/or inertial position information may be significant and 
cannot be neglected [Ref. 9]. In addition, the guidance and 
control laws must be as fast as possible in order to ensure 
accurate path keeping characteristics. The navigational 
position information time lags, as well as the dynamic lags 
due to the vehicle inertia, however, set a limit on the 
vehicle reaction time. Therefore, stability of the combined 


scheme becomes an issue which requires analysis. 


B. OBJECTIVE 

Previous studies have established boundaries for guidance 
and control laws in the horizontal [Ref. 10] and vertical 
plane [Ref. 12] maneuvering along straight line paths, as well 
as curved segments [Ref. 11]. The most important assumption 
in those results was the existence of instantaneous positional 
information updates when needed. In this work we relax the 
latter assumption. Stability boundaries are computed 
parametrized by the amount of positional information time lag. 
Results are presented based on eigenvalue and frequency 
response techniques. This thesis builds on the linear 
Stability analysis performed in [Ref. 13] which can predict 


the shape of the stability boundaries. In this work, vehicle 


motions in the vicinity of the corresponding boundary are 
assessed using nonlinear bifurcation theory [Ref. 5]. It is 
shown that the loss of stability occurs as generic Hopf 
bifurcations, where upon loss of stability of straight line 
motion a family of periodic solutions appears. Nonlinear 
expansions utilizing center manifold approximations and 
integral averaging techniques, reveal that these Hopf 
bifurcations can occur either in supercritical or subcritical 
forms. In the supercritical case, a stable family of limit 
cycles is generated immediately after the nominal motion 
becomes unstable. In the subcritical case, however, the 
resulting limit cycles are initially unstable and they are 
generated even before the nominal motion loses its stability. 
In the latter case, the domain of convergence of straight line 
motion becomes increasingly smaller as the stability boundary 
is reached. As a result, the methods developed in this work 
characterize the degree of confidence of the computed 
Stability boundaries by examining stability under finite 
external disturbances. A first order memory scheme, which can 
be easily implemented in real time, is suggested in order to 
expand the region of stability of straight line motion. All 
computations are performed for the NPS autonomous underwater 
vehicle for which a set of geometric properties and slow 
motion hydrodynamic derivatives is available [Ref. 1]. Unless 
Otherwise mentioned, all results are presented in standard 


nondimensional form. Equation development presented in 


[Ref.13] has been used in chapter II, and with modifications 


in chapter IV. 


C. THESIS OUTLINE 

Chapter II presents the formulation of vehicle equations 
of motion and the criterion for determining the region of 
guidance/control stability in the presence of a position 
information time lag. 

Chapter III examines the vehicle’s control stability 
through the use of a Hopf bifurcation analysis. 

Chapter IV examines the effect on the region of stability 
determination by the use of a memory model which incorporates 


the two previous vehicle positions. 


II. PROBLEM FORMULATION 


A. EQUATIONS OF MOTION 

The mathematical model which describes vehicle motion in 
the horizontal plane consists of nonlinear sway and yaw 
equations of motion. A moving coordinate frame which is fixed 
at the vehicle’s geometric center gives Newtonian equations 


of motion which are 


TV SUT aoe Vr“) = V (Za) 


I,r+MX,(Vtur) -my,vr=N (ec 


where u is constant forward speed; v and r are the relative 
Sway and yaw velocities of the moving vehicle relative to the 
water; m is the mass of the moving vehicle; Xg, yg are the 
respective lateral and longitudinal positions of the center of 
Gravity; and Y, N represent the total excitation sway force 
and yaw moment, respectively. The vehicle’s added mass, 
damping, and viscous drag due to motion through the water are 
represented by quadratic drag terms and memoryless polynomials 
inv and r. Y and N can be represented as the sum of these 


terms by 


Y= se Y,r+ eee v+Y_ur) + BiG, uv 


- Of" c,h D Blak Sai cat 2 2yyu2t 


tail “Y [v+Er| 


ae TW, r+ eu (N,v+N_ur) ee, uv 


er oe eg ivaea le 3n7 772 
£fr ic pf [veer] V+Er| Edt+f A N,u 6 


where 1 is the vehicle length, Y N, represent (parser 


a’ 


derivatives of Y and N with respect to a, C,. 1S (EHemieme 


y 
coefficient, and 6 is the rudder angle. 

At relatively high speeds and low angles of attack the 
steering response is predominantly linear. In low speed 
maneuvering or hovering conditions the cross flow integral 


drag term becomes significant. 


Vehicle yaw rate is expressed by 


wH=r (2 33)) 


and inertial position rate is expressed by 


y=usinyt+vcosyp (24) 


where y is the heading angle and shown in Figure 1. 








Figure 1. Vehicle Geometry and Definition of Symbols. 


Taking equations (2.1), (2.2), and (2.3) as a set cf tare 
coupled linear differential equations, and a linearized form 
of equation (2.4), the final set of equations of motion for 


steering control are 


7 (2352) 
V=a,,uvta,,urt+b,u*s (2.5b) 
f=a,,uv+a,,ur+b,u*s (2251) 

yrupty ( 22a) 


at any nominal forward speed. 


B. GUIDANCE AND CONTROL 
1. Guidance Scheme 

An orientation based command scheme is used where the 
commanded vehicle heading angle y. 18 a function of theme 
of sight angle o between the actual vehicle position and a 
reference position on the nominal path located at a constant 
lookahead distance d. Schematically, this is presented in 
Pagure =i. 


The line of sight angle o is defined by 


Sas 2.6 
tano - ( ) 


The autopilot will be called upon to orient the 
meneete tO the commanded heading angle y,, which, in pursuit 
Guidance, will equal the line of sight angle. This defines y, 


as 


Wo=-arctan+ (227) 


For relatively small angles 


eae (7 ec) 


2. Controller Design 
The steering control governing equations can be 


represented in the form 


x =Ax+béd (23) 


where the state vector equation is 


ea Weve (2.10) 


Linear full state feedback is introduced in the form 


\O 


S=k, (W-.) +k v+k,z (2m 


The gains k,, kj, and k, are computed by pole placement 
to give the desired dynamics to the closed loop system. The 
vehicle’s longitudinal axis is pointed toward the desired 
heading by the difference (¥-y,) in the control law. 

With a dimensionless time constant t, the general form 


of the characteristic equation is 


AP +a, A%+a,A+0,=0 (2qe2e 
where 
275 eS alk 
Be 


The controller gains are computed from 


a . 
ee (2. 139) 
*  (b,a,,-b,a,,) uw? 


(b,a,,-b,a,,) u*k,+ (b,a,,-b,a,,) u*k, 


(2 sali) 
ee pera c= (el each ela.) el 


but ko +b, 7k = =05 (eee ae (2 13@) 


160 


These gain values are continuously updated due to changing 
nominal forward speeds. Figure 2 shows the three controller 


gains at unit forward speed versus the system time constant. 


C. TIME LAG 

All of the required state information for vehicle control 
is available at sufficient rates with the possible exception 
of the y position information. Due to time requirements for 
Peeuwetion of sonar returns and inertial navigation 
information, there may be a time lag in the y position 
information. 

This time lag T (sec) is introduced to the guidance 


control law and the commanded rudder angle y¥, becomes 


p.=-arctanZte 7 (2.14) 
For small angles 


Te ea (2.15) 


and the resulting linearized control law becomes 


b=k,yrk, LED sk vekgr (utc 


Le 


CONTROLLER GAINS 





te 


Figure 2. Controller Gains versus t, at unit forward speed. 


(k,=solid, k.=dashed, k,=dotted) 


12 


The linearized equations of motion become 


W=r 


Via OU onl rel ita. t1-K,) Vv: 


k 
+ (a,,u+b,u*k,) r+b,ut—y( fat) 


© = Deck Wri(a,,utD,ucK,) v 


k 
ae (eh aweniulas) r+b,u*—-y( a2) 


yeupty 


D. STABILITY ANALYSIS 


(25 Fas) 


(el 7D) 


(2a 7e) 


(2219 dy) 


Control law stability is determined by the eigenvalues of 


the matrix A from the linear system 


xX=AX 


with the state vector 


= Vt py 


Le 


(28) 


(2.19) 


where the state equations have been linearized around straight 


Line, motion. The characteristic equation of (2.18) is a 


quartic of the form 


A*+BA?+CA*+DA+E=0 (2.20) 


where the coefficients B, C, D, and E are functions of vehicle 
properties, controller gains, and lookahead distance. 
The characteristic equation will give a pair of complex 


conjugate roots which cross the imaginary axis under the 


following conditions. 


BCD-B*E-D*=0 (2522a9) 


bis 


—=>0 (22225) 


Minimum lookahead distance d_¥;, for Stability (eam 


determined by translating these conditions to 


a,d-tanaay— 0 (2.22a) 


(222215) 


where 


14 


1 Oe aes eZ) 


(2 223) 


Taube 


a = ab P1G227 241272) [by * (by) a22~ B24; 7b,) ul @, (21236) 
: (b,a,,—b,a2,) *u 


Analysis of the system when operating with an information 
time lag requires approximation by either a first order Taylor 


expansion 
Wed) = y— Dy. (2.24) 
a second order Taylor expansion 


2. 
y(t-T) #y-Ty+ Geos 


a third order Taylor expansion 


2 3 
y(t-T) 2y-Ty+ 9 (Ze26) 


Or a frequency response analysis using the Nyquist criterion. 


Figure 3 shows the region of stability given by each of these 


i95) 


methods with a one second (dimentionless) time lag. It can be 
seen that the agreement is, in general, very good among the 
three Taylor series approximations, especially as the time 
COnetant £2) 1S sneredsecr The third order approximation 
coincides with the exact solution from the frequency response 
(Nyquist criterion) method. The first order method requires 
the highest value of d for stability (at ag iven sce ee 
therefore is the most conservative method for design use. For 
this reason the first order method is used in the next chapter 
to study the dynamics of the system as the critical stability 


boundary is crossed. 
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Figure 3. 
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ee 


Taylor expansion and Nyquist stability analysis 


with a one second time lag. (Taylor first 


order=solid, Taylor second order=dashed, Taylor 


third order and Nyquist=dotted) 


Ay 


III. HOPF BIFURCATION 


A. INTRODUCTION 

It can be numerically confirmed that in all cases of 
Stability loss of the previous chapter, one pair of complex 
conjugate eigenvalues of the corresponding eigenvalue problem 
crosses transversely the imaginary axis. A situation like 
this in which a certain parameter is varied such that the real 
part of one pair of complex conjugate eigenvalues of the 
linearized system matrix crosses zero, will result in the 
system leaving its steady state in an oscillatory manner. 
This loss of stability is called Hopf bifurcation and 
generically occurs in one of two ways, supercritical or 
Subcritical. In the supercritical case, stable limit cycles 
are generated after the nominal straight line motion loses its 
Stability. The amplitudes of these limit cycles are 
continuously increasing as the parameter distance from its 
critical value is increased. For small values of this 
criticality distance the resulting limit cycle is of small 
amplitude and differs little from the initial nominal state. 
In the subcritical case, however, unstable limit cycles are 
generated before the nominal state loses its stability. 
Therefore, depending on the initial conditions it is possible 


to diverge away from the nominal straight line path and 


die) 


converge toward a different steady state even before the 
nominal motion loses its stability. In many cases this new 
steady state is another limit cycle of considerably higher 
amplitude. This means that in the subcritical Hopf 
bifurcation case the domain of attraction of the nominal state 
is decreasing and in fact it shrinks to zero as the critical 
point is approached. Random external disturbances of 
sufficient magnitude can throw the vehicle off to an 
oscillatory steady state even though the nominal state may 
still remain stable. After the nominal state becomes 
unstable, a discontinuous increase in the magnitude of motions 
is observed as there exists no simple stable nearby attractors 
for the vehicle to converge to. Distinction between these two 
qualitatively different types of bifurcation is, therefore, 
essential in design. The computational procedure requires 
higher order approximations in the equations of motion and is 


the subject of the next section. 


B. COMPUTATIONS 


The vehicle equations of motion are written in the form 


we=r (3 .1a) 


V=a,,uvta,,ur+b,u7d (Seb) 


ies 


Y=a,,uVta,,ur+ pou 6 (3 ey 


y=usinyt+vcosyp (321a) 


In order to capture the effect of rudder saturation we 


write the rudder angle 6 in the form 


5 
sat tanh (=> 


6 =6 ) (320 





sat 


where 6.., 1s the saturation limit of 6, typically Seteqieamee 
radians. Equation (3.2) is used instead of a hard saturation 
function because of its smoothness, which is important for the 


computation in this section. 6, is the linear controia wae. 


the form 
6, =k, (W-.) +k v+tk yr (3.33 
where 
SEV A GaT) eae LY, 3.4 
tany. = 5 ( ) 


using the first order approximation for y(t-T). 


The state equations (3.1) through (3.4) are written in the 


GOMpacce. Eom 


ZU 


x= f(x) (2 3) 


where f(x) is a nonlinear function of the state variables 


vector 


x= (vr, y) (35 6)) 


Expanding f(x) in Taylor series, we can write (3.5) in the 


form 


X= AXx+g (xX) (a) 


where A is the linear system matrix and g(x) contains all the 
leading nonlinear terms of f(x). Due to the port/starboard 
symmetry of our problem, g(x) contains only third order terms. 
Defining T as the matrix of eigenvalues of A and applying the 


transformation 


the state equation is transformed into its canonical form 


Z-T tATzZ+1. 9g ( Tz) (329) 


At the Hopf bifurcation point 


fa Al, 


0 -w, 0 0 


riar|°o ° ° © (33a 
0 OSS 0 
0 8) 
where p and q are negative and wy, is positive. In the new 
system of coordinates 
2-1 (3 opr) 


the dynamics of the system are generated by a reduced two 
dimensional system z, and z5, since the coordinates Zi manam7, 
correspond to the eigenvalues p and g and are asymptotically 
stable. These stable variables zZ,, 2, Can be expres eegga ms 
function of the critical variables z,, zZ5; these funeewena 
expressions are at least of third order. Therefore, the 
stable coordinates Z,, 2, do not influence the thizd@owge. 
expansion of (3.7). Using the above expression, we can write 


the two dimensional system in Z,, Z,) as 


—o 3 2 iz 3) 


ae 3 2 2 3 
Z, = O24 725, 21 she ee ee (3. TZ 


where the coefficients sae Ee uelS complicated expressions that 


are derived from (3.9) 


Ze 


Equations (3.12) are valid at exactly the Hopf 
bifurcation point. For values of the parameter d close to the 


bifurcation point, they are written as 


pene nw fz +2, 2.) (aie 3 a) 


Pee cunt zee r (2. 2.) (Serb) 


‘ 


where a’, w’ are the derivatives of the real and imaginary 
parts of the critical pair of eigenvalues evaluated at the 
bifurcation point; € is the difference of the parameter d from 
@eemeritical value; and the nonlinear function F,, F, remain 


the same as in (3.12), 


= 5 2 2 3 


ol a 2 2 3 
15, aaa) Gn ol ee eee ee AP eee (3.14b) 


If we introduce polar coordinates in the form 


Zia COs¢ (3154) 


Z,=RsinO (37S) 


equations (3.13) become 


Ze 


R=a’eR+F, (R,8) cos0+F, (R,8) sind (3.16a) 


RO = (w,+w’e) R+F, (R, 8) cos6-F, (R, 6) sind (3/16) 


Equation (3.16a) yields 


R=a’eR+P(8) R? (30) 


where P(@+27)=P(6). If (3.17) is averaged over one cycle in 


6, we get an equation with constant coefficients 


R=«a’eR+KR? (37 ee 


where 


K= = [°"P(8) ® (3.19) 
0 


Carrying out the indicated integration in (3.19) we obtain 


K=>— (3277. ett toe Gere 2101) 


cole 


Existence and stability of limit cycles can be determined 
by analyzing the equilibrium points of the averaged equation 
(3.18), which correspond to periodic solutions in 2), §a@5ee 
can be seen from (3.15). If we examine equation (3.18) we can 


see that: 


24 


di) If a’>0, then: 


(a) If K>0 then unstable limit cycles coexist with the 
stable equilibrium for é<0. 


(b) If K<0 then stable limit cycles coexist with the 
unstable equilibrium for ¢é>0. 


Zoe Lf w’<0, then: 


(a) If K>O then unstable limit cycles coexist with the 
Stagle equilibrium for é>0. 


(b) If K<O then stable limit cycles coexist with the 
unstable equilibrium for é<0. 

Therefore, we can see that computation of the above 

nonlinear coefficient K can distinguish between the two 


different types of Hopf bifurcation: 


@® Supercritical if K<0 
Beoubcritical if K>0 


(Ref. 10). 


C. RESULTS 

Values of the coefficient K have been calculated for 
several different values of a y position information time lag, 
Over a span of system time constants, using the Fortran 
program presented in appendix A. This program has also been 
used to determine the vehicle period of oscillation in the 
Same conditions. Figure 4 shows the behavior of the K 
coefficient at three different time lags, and Figure 5 shows 


periods of oscillation at the same three time lags. It can be 


ZS 


Figure 4. 





te 


Coefficient K versus t, at three values of time lag 


Te: (0.5 sec=solid, 1.0 sec=dashed, 1.5 sec=dotted) 
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PERIOD 





tc 


Figure 5. Vehicle oscillation period versus t, for three 
values of time lag T. (0.5 sec=solid, 1.0 


sec=dashed, 1.5 sec=dotted) 
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seen that supercritical bifurcations are ensured for 
sufficiently high values of the control time constant Comes 
t, less than a certain critical value, the corresponding@ieps 
bifurcations are subcritical. This situation should be 


avoided in practice. 


D. SIMULATIONS 

The vehicle path has been simulated by using the Euler 
integration method coded in a Fortran program presented in 
appendix B. The vehicle control law (2.17) as well as the 
control gains (2.13) are used. The program has been run with 
input values for system time constant, y position information 
time lag, and lookahead distance. The vehicle was given an 
initial nominal y offset of half a shiplength and a nominal 
forward speed. The resulting plots of y position against time 
show the vehicle’s oscillatory path either converge to the 
commanded path or diverge. 

The results of these simulation runs were compared to the 
results of the Hopf bifurcation computations in two ways. The 
period of oscillation observed in the vehicle path was 
compared to the period predicted by the Hopf program. 
Additionally, the vehicle stability, determined by convergence 
to the commanded path, was compared to the K coefficient sign 
prediction of stability which was determined by the Hopf 
computations. Figure 6 shows a vehicle path in stable 


conditions. Figure 7 shows a vehicle path in unstable 


Zo 


Genditions. These simulations were chosen in the 
Supercritical bifurcation case and it can be seen that the 
period of oscillation observed from Figures 7 agrees very well 


with the theoretical results of Figure 5. 


Zo 


Figure 6. 


20 30 40 50 60 70 80 90 100 


Vehicle path in stable conditions where time 
constant t, = 1.0 sec, time lag T = 1.0 sec, and 


lookahead distance d = 2 shiplengths. 
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Figure 7. 


20 30 40 00 60 70 80 90 100 


Vehicle path in unstable conditions where time 
constant t, = 1.0 sec, time lag T = 1.0 sec, and 


lookahead distance d = 1.25 shiplength. 
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EV; STABILITY ENHANCEMENT 


A. TWO-POINT FORMULA 

The time lag of y position information discussed in 
chapter II has thus far been represented by a single piece of 
information used in the control law with some delay assigned 
to it. In an attempt to improve stability with regard to this 
delay, the use of the previous two positions in the control 
was examined. 

Equation (2.15) shows the representation of the time 
delay, and equation (2.16) shows its effect on the control 
law. A straight line approximation is applied to the previous 
two poSitions as shown in Figure 8. The general y position is 


defined by 


Ve Giese yeas (4.1) 


The two previous positions become 


yY, =a, 0, +a, (4.2a) 


Y2,=a,C, +a, (422) 


The coefficients a, and a, are stated as 
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Figure 8. Straight line approximation using two previous 


positions. 


ae 


a, 2 
C,~t, 
,= eee: (4 
t,-C, 
The previous time values are represented as 
CG, = been (4 
te 6 (4 


.3a) 


ycile) 


.4a) 


.4b) 


where T is the time lag associated with the y position 


information. The previous positions, in terms of the previous 


time values are 


Ve i ems 1) (4. 


Vo ae) (4. 


5a) 


= ))8)) 


Substitution and algebra gives a general position 


expression of 


y(t) =2y(t-T) -y(t-2T) (4.6) 
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B. STABILITY ANALYSIS 


The Laplace transform of equation (4.6) is 


Wace) =v (Bes =—e 727") (a7) 


The control law becomes 


k 
SakiWekjvek,r+—y(2e"-e 7%) (4.8) 


and the linearized equations of motion become 


WHr (4.9a) 


Ve etoile det uk) V+ (a,,u+tb,u-k,) rz 


k (4.9b) 
+b, u*—y (2e -e 27) 
QS a att. ok | Vr (aout bu kr , ? 
49°C 
+b,u*@—y(2e M-e 27) 
y=uyt+v (45 oq) 
These motion equations reduce to the state space form 
po ALS (A 10) 
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The matrix form becomes 


0 0 iL 0 y 


k 
2 2 2 oye tet SAN a BUS 
1U°K, @),U+D UK, 2),U+ 2) Ue eee —q fe e ) 


k 
prec, ei Deslon Vicl a a,,U+ pou ie bu, (2e-Ts-e-2Ts) |r 
u al 0 0 





The characteristic equation of the form [A-Is]=0 is 


s*‘+A,S° +A, S°tA, Ss? (6, S "2b ste o (2c sec (aes 


where 


A, =~ (4),+42,) U- (bk, +b,k,) u* 


A, = (411422742421) U** (4,5, -42,b,) u*k, + (a,b, -a,,b,) u°k, —b,u*k, 


A, => (45,0) 32), 0) nue 


B,==b,u*k,=(4,s bose ee 


S15 


15h (elie eben) baangy 


The characteristic equation of the system is written as 


1+KG(s) =0 (4.12) 


where 


(BN +B Seenice 2-6 275) 


G(s) = = : (4.13) 
S(s°+A,S°+A,S+A,) 
is the transfer function, and we have denoted 
Ke (4.14) 
With s=jw, the phase angle is represented by 
o=ZG(jw) (42.175) 


where 


b= 2 (2e°9%-e@ 27) -/ (jw) +Z (-B,w?+B,jw+B,) 
SW AGW? Aegon Ay) 


a7 


1 =2sin(oT)+sin(207) _ ft ie 


o@=tan7 —t+tan™ 
2cos (wT) -cos (2WT) 2 B.-B.w? 
-w?-A,w a al 
-tan™? (eee 
-A,W*°+A, 


The Nyquist stability criterion states that at the solution of 


o(w) =-nt (4 Flee 


where w is the phase cross over frequency o,, the gain) mangan 


must equal 1 


|KG(jw,)|=1 (4.18) 


where K = D = 1/d. Equation (4.17) becomes 


“2511 (@)7) +Sini(Zong Bw 
tan ?}———___4 a - B etan t (++) 
2608S (,7) -Cos(2Oyn mee B,- Bow: 
W)+A,W 
=tan) (ae en 
A,WitA, 
(4.19) 
Rearranged as 
G) @)+A,W 
tan” (—+—+_ ) -tan? (++) = 
_T _ pan 1 28in (0,7) +8in (20,7) 
2 2cos (w,7) -cos(20,7) 


3g 


Taking the tangent of both sides and using the identity 


tan (x) -tan(y) 


eon eae) Can (y) 


then algebra and rearranging gives 
B,W, (A,-A,@7) ~ (By-B, 03) (-01+A,0,) _ 
(B,-B, 03) (-A,w}+A,) +B, W, (-w}+A,@,) (242:0) 


ZEOS(@,7-Cos (2), T) 
eeemn (ayy s it (26,7) 


The stability computations require an initial value for 
@,- This is accomplished by setting the time lag T=0, and 


using 


(B,-B,W{) (-A,W{+A,) +B,W, (-w{+A,0,) =0 (4.21) 


rearranged as 


(B,A,-B,) 01+ (B,A,-B,A,-B,A;,) ©) +B A, =0 (am22)) 


Setting the gain margin equal to 1 gives 


i! | (-B,0i+B, jw, +B) (2e 7-2 777) | 
me a ee 
e W,|-J01-A,01+A, 7, +A, | 


and d can be solved for by 
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B,-B,W*) +B*w, /5-4C08(0,T 
Ge ( 97 Po 1) 1 a (Anoa5 


2 A 
wy) (A,-A,@{) “+ (A,@,-@}) 


In order to facilitate computation, terms are grouped as 


follows 


a, =B,A,-B, 


a, =B,A,-B,A,- BA, 


Te SE 


Bb, =8,+8,A,—-B,A, 


B, =B,A,-B,A, 


Equation (4.20) becomes 


B,w°+B,0°+B,0 | 2cos (wT) -cos (2WT) 
BE 8 SSS (4.24) 
%,W4+a,W*+a, -2sin(wT) +sin(2WT) 


Put in the form 
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f(w) =0 


equation (4.24) becomes 


eyo +6.0°+6,0) (-2sin(wT) +sin(2T) )) 


ioe G) +e. (2cos (1) —-cos (207) ) ] =0 


MNeweon’s iteration is used as 


£(W,.;) 


@,=W® SS, 
kK k=1 y 
Gy) 


Where the function of w is 


f(w) = (PB, @°+B,0°+B,w) (-2sin(wT) ) 
+ (B,w°+B,w?+B,0) (sin(2WT) ) 
aG7@ +t0,0°+a.) (2cos (MT) } 
+ (a,0*+a,0%+a,) (cos (2WT) ) 


ama the derivative function is 


Ba (O) = (Sams 3 6 a@°+B,) (-25in (wT) ) 
=(p7O-4070 *B,@) (27TCos( a7) 
(Spor +3Rrotep.) (sim(2wT) ) 
+(B,w°+B,w°+B,0) (2Tcos(2WT) ) 
+ (4@,0°+2a,W) (-2cos (WT) ) 
+(a,w*ta,w*+a,) (2TSsin(WT) ) 

+ (4a,wW°+2a,W) (cos (2WT) ) 
-(a,0*ta,w*%+a,) (2Tsin(2wT) ) 
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(42'S ) 


(4.26) 


(4227) 


(4.28) 


(4.29) 


Appendix C presents the Fortran program used to determine 
stability characteristics as a function of system time 


constant and time lag. 


C. RESULTS AND DISCUSSION 

Values for the minimum lookahead distance required for 
controller stability have been computed as a function of the 
system time constant and the y position information time lag 
using the two-point formula. Figures 9, 10, and 11 show 
critical lookahead distance plotted against time constant at 
three values of time lag for both the two-point method as well 
as the single point (Nyquist) method discussed in chapter II. 
These figures show an overall decrease in required lookahead 
distance as a result of the use of the two-point memory model. 
An exception to this occurs at low values of t.. Howe vce 
these values are to be avoided in practice, as explained in 
the previous chapter. 

Figure 12 through 15 show the critical lookahead distance 
plotted against time lag for different values of the time 
constant uSing both the two-point and single point methods. 
Again the two-point method produces superior results to the 
Single point method, with the exception of small time 
constants and large time lags. The same information is 
summarized in Figure 16 through 19, where we show the ratio 
d,/d, (critical lookahead distance using the two-point method 


divided by the critical lookahead distance using the single 
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point method) verses time constant and time lag. It can be 
seen that the use of the two-point method can result in 


improvement in the region of stability of over 20%. 
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1.6 
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ee 


Figure 9. 


1 f.2 1.4 Bas) 1.8 Zz eee 


Critical lookahead distance d versus time constant 
t, with time lag T = 0.5 sec (two-point 


method=solid, single point=dashed). 
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Figure 10. Critical lookahead distance d versus time constant 


t, with time lag T = 1.0 sec (two-point 


method=solid, single point=dashed). 
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tc 





Figure 11. Critical lookahead distance d versus time constant 


t, with time lag T = 1.5 sec (two-point 


method=solid, single point=dashed). 
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1.8 


1.6 





Figure 12. Critical lookahead distance d versus time lag T 
with time constant t. = 0.8 sec (two-point 


method=solid, single point= dashed). 
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Figure 13. Critical lookahead distance d versus time lag T 


with time constant t. = 1.0 sec (two-point 


method=solid, single point= dashed). 
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3 Sere 


Figure 14. Critical lookahead distance d versus time lag T 
with time constant t. = 1.5 sec (two-point 


method=solid, single point=dashed) . 
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See 3.4 3.6 3.8 4 Ane Al 


Figure 15. Critical lookahead distance d versus time lag T 
with time constant t. = 2.5 sec (two-point 


method=solid, single point=dashed). 
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Figure 16. Ratio d./d, versus time constant t, with time 


lag = 1.0 sec. 
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d2/dl 


Figure 17. Ratio d./d, versus time constant t, with time 


lag = 1.5 sec. 
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d2/dl 


figure 18. Ratio d./d, versus time lag T with time constant = 


1.0 sec 
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d2/dl 


Figure 19. Ratio d,/d, versus time lag T with time 


constant = 1.5 sec. 
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Vv. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

Nonlinear bifurcation theory has been applied to the NPS 
autonomous underwater vehicle through a Hopf bifurcation 
analysis of the vehicle’s navigation guidance and control 
scheme. It has been shown that under normal conditions the 
vehicle will operate in the supercritical region, and 
therefore the unpredictable behavior associated with the 
subcritical case should not to be expected. 

Secondly, the control scheme was modified by the use of a 
Single stage memory model which incorporates the two previous 
vehicle positions in the linearized control law. It has been 
shown that the use of this memory model reduces the vehicle’s 
required lookhead distance for control stability within the 


normal operating range. 


B. RECOMMENDATIONS 

In the event that not all of the state information 
required is directly available through measurements, the 
control/guidance scheme would include an observer. The effect 
of the observer dynamics on the higher order nonlinear terms, 
and therefore on the bifurcation limit cycle stability, should 


be examined. 
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APPENDIX A 


PROGRAM HOPF.FOR 
Hopf Bifurcation Analysis 


Critical Point: Exact 
Third Order Expansions: First Order Approximation 


IMPLICIT DOUBLE PRECISION (A-H,O-Z) 

DOUBLE PRECISION K1,K2,K3,K4,L,NR,NV,NDRS,NDRB,IZ,MASS, 
NRDOT,NVDOT,K1P,K2P 

DOUBLE PRECISION M11,M12,M13,M14,M21,M22,M23,M24, 
M31,M32,M33,M34,M41,M42,M43,M44, 
NIL, N12, N13 NIA NZ Nee Nes Neo 
N31,N32,N33,N34,N41,N42,N43,N44, 
L21,L22,623,024 3 lo ooo, 
Lal) b42, Las La 


DIMENSION A(4,4),1T( 4,4), TD INV(4, 4) bVIC4)) hy ad 24 eee 
DIMENSION WR(4),WI(4),TSAVE(4,4),TLUD(4,4),IVLUD(4),SVLUD(4) 
DIMENSION ASAVE(4,4) 

OPEN (11,FILE=’HOPF.DAT’ ,STATUS='NEW’ ) 

OPEN (12,FILE=’ PERI.DAT’,STATUS=’ NEW’ ) 

OPEN (13,FILE=’ALPH10.DAT’ ,STATUS='NEW’ ) 

Vehicle Parameters 


WEIGHT=435.0 


IZ =45.0 

L =7.3 

RHO =1.94 

G =32.2 

XG =0.0104 
MASS =WEIGHT/G 
U =5.0 
RATIO =0.0 

DO =0.4 


WRITE (*,1006) 
READ (*,*) DO 


YRDOT =-0.00178*0.5*RHO*L**4 
YVDOT =-0.03430*0.5*RHO*L**3 
YR =+0.01187*0.5*RHO*L**3 
YV =-0.03896*0.5*RHO*L**2 
YDRS =+0.02345*0.5*RHO*L**2 
YDRB =+0.02345*0.5*RHO*L**2 
NRDOT =-0.00047*0.5*RHO*L**5 
NVDOT =-0.00178*0.5*RHO*L**4 
NR =-0.01022*0.5*RHO*L**4 
NV =-0.00769*0.5*RHO*L**3 
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IRO Sr-U.LSSIXNVU.VULZISGO*VU.O* AAU L"™"S 
ORB =+0.283*0.02345*0.5*RHO*L**3 


4 =(IZ—-NRDOT) *(MASS-YVDOT) - 

| (MASS *XG-YRDOT) * (MASS*XG-NVDOT) 
A11=((IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DH 
\12=( (1Z-NRDOT) *(-MASS+YR)- 

| (MASS*XG-YRDOT ) * (-MASS*XG+NR) ) /DH 
421=((MASS-YVDOT ) *NV-(MASS*XG-NVDOT) *YV) /DH 
A\22=( (MASS-YVDOT ) * (-MASS*XG+NR) - 
(MASS*XG-NVDOT) *(-MASS+YR) ) /DH 

311=( (1Z-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS) /DH 
312=((1IZ-NRDOT) * YDRB-(MASS*XG-YRDOT) *NDRB) /DH 
321=( (MASS-YVDOT ) *NDRS-(MASS*XG-NVDOT) *YDRS) /DH 


322=( (MASS-YVDOT ) *NDRB-(MASS*XG-NVDOT) *YDRB) /DH 


31=BB11+RATIO*BB12 
32=BB21+RATIO*BB22 


RITE (*,1005) 
SAD (*,*) Ti 
MITE (*,1001) 
SAD (*,*) Te 


»S=1.D-10 
"MAX=1000 


(L=0.5 
x0 50 M=1,100 
.D=TL*L/U 


ac=0.5 

» 40 K=1,100 

~—TCD=TC*L/U 

POLEC=1.0/TCD 

AD1=3.0*POLEC 

AD2=3.0*POLEC**2 

AD3=POLEC* * 3 

A1l=BB1*U*U 

B1=BB2*U*U 

C1=-AD1-(AA11+AA22)*U 
A2=(BB1*AA22-BB2*AA12 ) *U* * 3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
K1=AD3/((BB2*AA11-BB1*AA21 ) *U**3) 
C2=AD2-(AA11*AA22-AA12*AA21 ) *U**24+BB2*U*U*K1 
K2=(C1*B2-C2*B1)/(A1*B2-A2*B1 ) 
K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 


Al= (AA11*BB2-AA21*BB1 ) *U*U*U*K1 


A2= (AA11*AA22-AA12*AA21 ) *U*U+(AA11*BB2-AA21*BB1 ) *U*U*U*K34+ 


(AA22*BB1-AA12*BB2) *U*U*U*K2—-BB2*U*U*K1 
A3=- (AA11+AA22) *U-(BB1*K2+BB2*K3 ) *U*U 
BO= (AA11*BB2-AA21*BB1 ) *U*U*U*U*K1 
Bl=—(AA12*BB2-AA22*BB1 ) *U*U*U*K1—-BB2*U*U*U*K1 
B2=-BB1*U*U*K1 


Compute Initial Approximation For Omega (TL=0) 


SW 


BOMB SIA SHB ns 
BQ=B1*A2-B2*A1-B0*A3 

CQ=BO*Al 

DET=BQ**2-4.D0*AQ*CQ 

IF (DET.LT.0.D0) GO TO 40 
ZT1=(-BQ+DSQRT(DET) )/(2.0D0*AQ) 
ZT2=(-BQ-DSQRT( DET) )/(2.0D0*AQ) 
IF (ZT1.GE.0.D0) OM=DSQRT(ZT1) 
IF (ZT2.GE.0.D0) OM=DSQRT(ZT2) 
ALPHA1=AQ 

ALPHA2=BQ 

ALPHA3=CQ 

BETA1#=-B2 

BETA2=B0+B2*A2-B1*A3 
BETA3=B1*A1-BO0*A2 


Compute Exact Omega Using Newton’s Iteration 


Cerise! Sari ae Ge Cae tay Coe) 


OMOLD=OM 

93 F=(BETA1*OMOLD* *5+BETA2* OMOLD* * 3+BETA3*OMOLD) *DSIN( OMOLD*TLD) 
+ (ALPHA1 *OMOLD** 4+ALPHA2 *OMOLD* * 2+ALPHA3 ) *DCOS ( OMOLD*TLD) 
FPRIME=(5.D0*BETA1 *OMOLD** 4+3.D0*BETA2*OMOLD**2+BETA3) 
*DSIN(OMOLD* TLD) +( BETA1 *OMOLD* *5+BETA2*OMOLD* * 3+BETA3*OMOLD) 
« TLD*DCOS ( OMOLD*TLD) +(4.D0*ALPHA1 *OMOLD* * 3+2.DO*ALPHA2*OMOLD ) 
*DCOS ( OMOLD* TLD) — (ALPHA1 *OMOLD* * 4+ALPHA2* OMOLD* * 2+ALPHA 3 ) 
x*TLD*DSIN(OMOLD*TLD) 

IF (FPRIME.EQ.0.D0) STOP 1112 

IF (F.EQ.0.D0) GO TO 92 

OMNEW=OMOLD-F/FPRIME 

OMDIF=DABS( OMNEW-OMOLD) 

IF (OMDIF.LE.EPS) GO TO: 92 

OMOLD=OMNEW 

IT=IT+1 

IF (IT.GT.ITMAX) STOP pili 

GO TO 93 


Q 


RQ Qt Mr M 


‘92 OM=OMNEW 
XDNUM=DSQRT( (BO-B2*OM**2) **2+(B1*OM) **2) 
XDDEN=OM*DSORT( (A1-A3*OM**2) **2+(A2*OM-OM** 3) **2) 
XD=XDNUM/XDDEN 


Start Hopf Bifurcation Analysis 


Evaluate Nonlinear Rudder Expansion Coefficients 


qQagaan 


K1P=K1-K1*TLD*U/XD 
K2P=K2-K1*TLD/XD 


DPPV=- (1.D0/(3.D0*D0**2))*3.D0*KIP*KIP*K2P 
& + 0.5D0*K1*TLD/XD + 3.D0*K1*(TLD*U) **3/(3.D0*XD**3) 
DPVV=- (1.D0/(3.D0*D0**2) )*3.D0*KIP*K2P*K2P 
& + 3.D0*U*TLD*TLD*TLD*K1/(3.D0*XD**3) 
DPPR=- (1.D0/(3.D0*D0**2) )*3.DO*KIP*K1P*K3 
DPRR=- (1.D0/(3.D0*DO0**2) )*3.DO*KIP*K3*K3 
DPPY=- (1.D0/(3.D0*D0**2) )*3.DO0*KIP*K1P*K1/XD 


& — 3. D0*TLD*TLD*U*U*K1/(3.D0*XD**3) 

DPY Y= (1.D0/(3.D0*DO**2) )*3.DO*KIP*K1*K1/(XD**2) 
& 4 3, D0*TLDwUeKI (2. D0= xo. ts) 

DVVR=- (1.D0/(3.D0*D0**2) )*3.D0*K2P*K2P*K3 

DVRR=- (1.D0/(3.D0*D0**2) )*3.D0=K2P*K3"K3 


53 


— Se 


DVYY=- 

+ 
DRRY=- 
DRYY=- 
DPVR=- 
DPVY=- 
DPRY=- 
DVRY=- 
DPPP=- 

+ 
DVVV=- 

2b 
DRRR=- 
DYYY=- 


a. DO*K1*TLD*TLD/(3. DO*XD**3) 
(1.D0/(3.D0*D0**2))*3.D0*K1*K1*K2P/(XD**2) 
3.DO*TLD*K1/(3.D0*XD**3) 

(1.D0/(3.D0*D0**2) )*3.D0*K1*K3*K3/XD 
(1.D0/(3.D0*D0**2) )*3.D0*K1*K1*K3/(XD**2 ) 
(1.D0/(3.D0*D0O**2))*6.D0*KIP*K2P*K3 
(1.D0/(3.D0*D0**2))*6.D0*KIP*K1*K2P/XD 

6 .DO*TLD*TLD*U*K1/(3.D0*XD**3) 
(1.D0/(3.D0*D0O**2) )*6.D0*KIP*K1*K3/XD 
(1.D0/(3.D0*D0O**2))*6.D0*K1*K2P*K3/XD 
(1.D0/(3.D0*DO**2) )*1.D0*KIP*K1P*K1P 
K1*TLD*U/(6.DO*XD) + (K1*(TLD*U)**3)/(3.D0*XD**3) 
(1.D0/(3.D0*D0**2))*1.D0*K2P*K2P*K2P 
K1*(TLD**3)/(3.D0*XD**3) 

(1.D0/(3.D0*D0O**2) )*1.D0*K3*K3*K3 
(1.D0/(3.D0*D0**2) )*(K1/XD)**3-K1/(3.D0*XD**3) 
K1/(3.D0*XD**3) 


Evaluate Transformation Matrix of Eigenvectors 


PSI =0.D0 
Y =0.D0 
V =0.D0 
DPSI=K1P 
DV =K2P 
DR =K3 
DY =K1*XD/(XD**2+Y**2) 
A(1,1)=0.0D0 
A(1,2)=0.0D0 
mA(1,3)=1.0D0 
(1,4)=0.0D0 
(2,1)= BB1*U*U*DPSI 
(2,2)=AA11*U+BB1*U*U*DV 
(2,3)=AA12*U+BB1*U*U*DR 
(2,4)= BB1*U*U*DY 
(3,1)= BB2*U*U*DPSI 
A(3,2)=AA21*U+BB2*U*U*DV 
A(3,3)=AA22*U+BB2*U*U*DR 
A(3,4)= BB2*U*U*DY 
A(4,1)=U*DCOS(PSI)-V*DSIN( PSI) 
A(4,2)= DCOS(PSI) 
A(4,3)=0.0D0 
A(4,4)=0.0D0 
pO il I= 


ASAVE(1.3)*A(I,3) 
CONTINUE 
CONTINUE 
BALL RG(4,4,A,WR,W1,1,22Z,1V1,FV1,IERR) 
CALL DSOMEG(IEV,WR,WI,OMEGA, CHECK) 


OMEGA0= 
=1,4 


moO 5 I 


OMEGA 


Mer 1)= 222(1,1EV) 
ete, 2)=-222Z(1,1EV+1) 
CONTINUE 
MF (IEV.EQ.1) GO TO 13 
IF (IEV.EQ.2) GO TO 14 
IF (IEV.EQ.3) GO TO 15 
STOP 3004 
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02 0 


aqgaaAQ 


16 
17 


QR QM BW M1 MS LK 


T(1, 3 eo 22(1,1) 

T(I,4)=Z222(1,4) 

CONTINUE 

GOr To. 17 

DO 7 I=1,4 
T(1,3)=Z22(1,1) 
T(1,4)=Z222Z(1,2) 

CONTINUE 

GO TO 17 

DO 16 I=1,4 
T(1,3)Z2Z22(1,3) 
T(I,4)=222(1,4) 

CONTINUE 

CONTINUE 


Normalization of the Critical Eigenvector 


CALL NORMAL(T) 


Definition of Mij 


Definition of Lij 


L21= DPPV*M11*M11*M21 
+ DPRR*M11*M31*M31 
+ DVVR*M21*M21*M31 
+ DVYY*M21*M41*M41 
+ DPVR*M11*M21*M31 
+ DVRY*M21*M41*M31 
+ DRRR*M31*M31*M31 
PPV = M11*M11*M22 + 2 
PVV = M12*M21*M21 + 2 
PPR = M11*M11*M32 + 2 
PRR = M12*M31*M31 + 2 
PPY = M11*M11*M42 + 2 
PYY = M41*M41*M12 + 2 
VVR = M21*M21*M32 + 2 
VRR = M22*M31*M31 + 2 
VVY = M21*M21*M42 + 2 
VYY = M22*M41*M41 + 2 
RRY = M31*M31*M42 + 2 


oe 
+ 
+ 
+ 
+ 
+ 
+ 


DPVV*M11*M21*M21 
DPPY*M11*M11*M41 
DVRR*M21*M31*M31 
DRRY*M31*M31*M41 
DPVY*M11*M21*M41 
DPPP*M11*M11*M11 
DYYY*M41*M41*M41 


-O*M11*M12*M21 
~O*M11*M21*M22 
~O*M11*M12*M31 
-O*M11*M31*M32 
~O*M11*M12*M41 
~O*M11*M41*M42 
-O*M31*M21*M22 
-O@M3T*ns24M2r 
~O*M41*M21*M22 
~0O*M41*M42*M21 
~0*M41*M31*M32 
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DPPR*M11*M11*M31 
DPYY*M11*M41*M41 
DVVY*M21*M21*M41 
DRYY*M31*M41*M41 
DPRY*M11*M41*M31 
DVVV*M21*M21*M21 


VRY = 


ree More Mol + 2LUeMSl e Ma 2*Msl 
M11*M21*M324+M31*(M11*M22+M12*M21 ) 
M11*M21*M42+M41*(M11*M224+M12*M211 ) 
M11*M41*M324+M31*(M11*M424+M12*M41) 
M21*M41*M324+M31* (M21*M424+M22*M41) 
3.0*M11*M11*M12 

3.0*M21*M21*M22 

3.0*M31*M31*M32 

3.0*M41*M41*M42 


L22=DPPV* PPV+DPVV* PVV+DPPR* PPR+DPRR* PRR+DPPY*PPY+DPYY*PYY 
+DVVR* VVR+DVRR*VRR+DVVY* VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR* PVR+DPVY*PVY+DPRY* PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 


+DRRR* RRR+DYYY*YYY 
PPV = M12*M12*M21 + 2.0*M11*M12*M22 
PVV = M11*M22*M22 + 2.0*M12*M21*M22 
PPR = M12*M12*M31 + 2.0*M11*M12*M32 
PRR = M11*M32*M32 + 2.0*M12*M31*M32 
PPY = M12*M12*M41 + 2.0*M11*M12*M42 
PYY = M11*M42*M42 + 2.0*M12*M41*M42 
VVR = M22*M22*M31 + 2.0*M21*M22*M32 
VRR = M21*M32*M32 + 2.0*M22*M31*M32 
VVY = M22*M22*M41 + 2.0*M21*M22*M42 
VYY = M21*M42*M42 + 2.0*M22*M41*M42 
RRY = M32*M32*M41 + 2.0*M31*M32*M42 
RYY = M31*M42*M42 + 2.0*M32*M41*M42 
PVR = M12*M22*M31 + M32*(M11*M22+M12*M21) 
PVY = M12*M22*M41 + M42*(M11*M22+M12*M21) 
PRY = M12*M42*M31 + M32*(M11*M42+M12*M41) 
VRY = M22*M42*M31 + M32*(M21*M424+M22*M41) 
PPP = 3.0*M11*M12*M12 
VVV = 3.0*M21*M22*M22 
RRR = 3.0*M31*M32*M32 
YYY = 3.0*M41*M42*M42 


L23=DPPV*PPV+DPVV* PVV+DPPR* PPR+DPRR* PRR+DPPY*PPY+DPYY*PYY 
+DVVR* VVR+DVRR*VRR+DVVY*VVY+DVYY* VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 
+DRRR*RRR+DYYY*YYY 


c 
to 
ras 


t++eeett i 


L31=L 
L32=L 
L33=L 
L34=L 


L21=L 
L22=L 
L23=L 
L24=L 
L31=L 
1L32=L 


DPPV*M12*M12*M22 
DPRR*M12*M32*M32 
DVVR*M22*M22*M32 
DVYY*M22*M42*M42 
DPVR*M12*M22*M32 
DVRY*M22*M42*M32 
DRRR*M32*M32*M32 


ZA 
Os 
Z3 
24 


21*BB1*U*U 
22*BB1*U*U 
23*BB1*U*U 
24*BB1*U*U 
31*BB2*U*U 
32*BB2*U*U 


DPVV*¥M12*M22*M22 
DPPY*M12*M12*M42 
DVRR*M22*M32*M32 
DRRY*M32*M32*M42 
DPVY*M12*M22*M42 
DPPP*M12*M12*M12 
DYYY*M42*M42*M42 
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DPPR*MI2*M1Z2*M52 
DPYY*M12*M42*M42 
DVVY*M22*M22*M42 
DRYY*M32*M42*M42 
DPRY*M12*M42*M32 
DVVV*¥M22*M22*M22 


aan 


“Fs Faery 


ry i FY 


€2¢02 02 


202 02 


N11=TINV( 


wIInwhy I Soe uey 
L34=L34*BB2*U*U 
L41=-0.5*M11*M11*(M214+U*M11/3.0) 


L42=-M11*(M12*M21+0.5*M11*M22+0.5*U*M12*M12) 
L43=-M12* (M11*M22+0.5*M12*M21+0.5*U*M11*M12) 


L44=-0.5*M12*M12*(M224+U*M12/3.0) 


Invert Transformation Matrix 


CONTINUE 
CONTINUE 
CALL DLUD(4,4,TSAVE,4,TLUD,IVLUD) 
DO 4 I=1,4 
LF (IVEUD(1).EO.0) Sitter 3003 
CONTINUE 
CALL DILU(4,4,TLUD,IVLUD,SVLUD) 
DO 8 I=1,4 
DO 9 J=1,4 
TINV(1, 09 )=TLUD(1, 0) 
CONTINUE 
CONTINUE 


Check Inv(T)*A*T 


IMULT=2 


IF (IMULT.EQ.1) CALL MULT(TINV,ASAVE,T) 


IF SCINULT BO. 0)  STOr 


Definition of Nij 


1, 
N21=TINV(2, 
N31=TINV(3, 
N41=TINV(4, 
N12=TINV(1, 
N22@TINV(2, 
N32=TINV(3, 
N42=TINV(4, 
N13=TINV(1, 
N23=TINV(2, 
N33=TINV(3, 
N43=TINV(4, 
N14=TINV(1, 
N24=TINV(2, 
N34=TINV(3, 
N44=TINV(4, 
R11=N12*L214N13*L314+N14*L41 
R12=N12*L224+N13*L32+N14*L42 
R13=N12*L234N13*L334+N14*L43 
R14=N12*L244N13*L344N14*L44 
R21=N22*L214+N23*L314+N24*L41 
R22=N22*L224N23*L324+N24*L42 
R23=N22*L234+N23*L334+N24*L43 
R24=N22*L244N23*L344N24*L44 
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Bvaluate ALlplla ana VUMega- 


DINC=0.001 
XDR =XD+DINC 
XDL #=XD-DINC 


XD =XDR 
PSI =0.D0 
Y =0.D0 
Vv =0.D0 
DPSI=K1P 
DV =K2P 
DR #=K3 


DY =K1*XD/(XD**2+Y**2) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

0.0D0 

BB1*U*U*DPSI 

AA11*U+BB1*U*U*DV 

AA12*U+BB1*U*U*DR 
BB1*U*U*DY 
BB2*U*U*DPSI 

21*U+BB2*U*U*DV 

22*U+BB2*U*U*DR 
BB2*U*U*DY 

U*DCOS(PSI)-V*DSIN(PSI) 

a 


& & & & WWW WddN DNDN 


= #8 = = = = & = = = = SS = 


LS 
1) 
2) 
3) 
4) 
1) 
2) 
3) 
4) 
1) 
2) 
3) 
4) 


Prop y> > yyy DD DD 
veuuéH "Be hun iu 


BALL RG(4,4,A,WR,W1,0,22Z,1V1,FV1,1ERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 
ALPHR=DEOS 


OMEGR=FREQ 
XD=XDL 

PSI =0.D0 
¥ =0.D0 
Vv =0.D0 
DPSI=K1P 
DV =K2P 
DR =K3 


DY =K1*XD/(XD**2+Y**2) 
A(1,1)=0.0D0 
A(1,2)=0.0D0 
)=1.0D0 
0.0D0 
BB1*U*U*DPSI 
AA11*U+BB1*U*U*DV 
AA12*U+BB1*U*U*DR 
BB1*U*U*DY 
BB2*U*U*DPSI 
AA21*U+BB2*U*U*DV 
AA22*U+BB2*U*U*DR 
BB2*U*U*DY 
U*DCOS(PSI)-V*DSIN(PSI) 
DCOS( PSI) 


A( 
( 
( 
( 
( 

A ( 

A( 

A( 
( 
( 
( 
( 
( .0D0 
( 


eae ee ee 


3 
4) 
,1) 
72) 
73) 
74) 
71) 
72) 
73) 
4) 
1) 
72) 
73) 
74) 


1S 
= 
0 
oO 
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(2 


erzervera 


CALL DSTABL(DEOS,WR,WI, FREQ) 
ALPHL=DEOS 
OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) /(XDR-XDL) 
DOMEGA=( OMEGR-OMEGL) /( XDR-XDL) 


Evaluation of Hopf Bifurcation Coefficients 


COEF1=3.0*R11+R13+R22+3.0*R24 
COEF2=3.0*R21+R23-R12-3.0*R14 
AMU2 =-COEF1/(8.0*DALPHA) 
BETA2=0.25*COEF1 
TAU2 =-(COEF2-DOMEGA*COEF1/DALPHA) /(8.0*OMEGAO ) 
PERD =2.0*3.1415927/OMEGA0 
PER =PERD*U/L 
X=XD/L 
K4=K1/XD 
WRITE (11,*) TG, COBEL 
WRITE (12,*) TC,PER 
WRITE (13,%*) TC,DALPHA 
TC=TC+0.015 
40 CONTINUE 
TL=TL+0.02 
50 CONTINUE 


1001 FORMAT 


(‘ ENTER TIME CONSTANT’ ) 
1005 FORMAT s( 
( 
( 


ENTER TIME LAG’ ) 

1006 FORMAT 

2002 FORMAT 
END 


SUBROUTINE DSOMEG(IJK,WR,WI,OMEGA, CHECK ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 
CHECK=-1.0D+25 
pO 1 I=1,4 
IF (WR(I).LT.CHECK) GO TO l 
CHECK=WR(I) 
IJ=I 
1 CONTINUE 
OMEGA=DABS(WI(IJ)) 
IF (WI(IJ).GT.0.D0) IJK=IJ 
TE 4W0 (1d) bTs 0. D0) TJKets-—1 
RETURN 


SUBROUTINE DSTABL(DEOS,WR,WI1,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI1(4) 
DEOS=-1 .0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO l 
DEOS=WR(1) 
IJ=I 
i CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS ( OMEGA) 
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ND 

,-UBROUTINE NORMAL (T) 

IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 
IMENSION 1T(4,4),TNOR(4, 4) 
FAC=T(1,1)**24+T(1,2)**2 

F (DABS(CFAC).LE.(1.D-10)) STOP 4001 
NOR(1,1)=1.D0 


mee ye T(1,1)*T(2,1)4+T(2,2)*T(1,2))/CFAC 
NOR(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
NOR(4,1)=(T(1,1)*T(4,1)4+T(4,2)*T(1,2))/CFAC 
‘NOR(1,2)=0.D0 
INOR(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
INNOR(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
NNOR(4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2))/CFAC 
© 1 I=1,4 

BO 2 J=1,2 

T(1I,J)=TNOR(I,J) 

CONTINUE 
“ONTINUE 
\ETURN 
‘ND 


—_ ee we wes we wes ws ws we ws ws ws ww ww ww we ww ww we ww we we ww we we ws ww ww ws ww we ww we ee i i es ss ee 
=fS2sS2=S2=S2S2S=S=—2S-=S-=S-2=->=2=fsSS=S-S_=S-_=S—25S°—Q°0°“[—_°—_S—S-S-S=S-_S-_SS-_S2SSS- SSS S—=S—SSSSS-S-_ SS SSS SSS2S SS. SS =] = 


UBROUTINE MULT(TINV,A,T) 

MPLICIT DOUBLE PRECISION (A-H,0-Z) 

YIMENSION TINV(4,4),A(4,4),T(4,4),A1(4,4),A2(4,4) 
O 1 I=1,4 


A(1I,K)*T(K,J)+A1(1I,J) 


CONTINUE 
‘ONTINUE 
m 6 I=1,4 

DO 7 J=1,4 

DO 8 K=1,4 
A2(1,J)=TINV(I,K)*A1(K,J)+A2(1I,J) 
CONTINUE 

CONTINUE 
‘ONTINUE 
ma ll I=1,4 

merre (*,101) (A(I,J),J3=1,4) 
[CONTINUE 
m 12 I=1,4 

Smite (*,101) (T(I,J),J=#1,4) 
SONTINUE 
0 10 I=1,4 

WRITE (*,101) (A2(1I,J),J=1,4) 
SONTINUE 
WRITE (*,101) A2(1,1) 
,ETURN 
"ORMAT (4E15.5) 
IND 


65 


& 


& 


& 


& 


APPENDIX B 


PROGRAM SIMU.FOR (VEHICLE PATH SIMULATION) 
REAL K1,K2,K3,L,NR,NV,NDRS,NDRB,1IZ,MASS, 
NRDOT,NVDOT,K1P,K2P 


OPEN (20,FILE=’SIM1.DAT’,STATUS=’NEW’ ) 
OPEN (30,FILE=’SIM2.DAT’,STATUS='NEW’ ) 


WEIGHT=435.0 


IZ =45.0 

L =] .3 

RHO =1.94 

G =#32.2 

XG =0.0104 

MASS =WEIGHT/G 

U =5.0 

RATIO =0.0 

DO =0.4 

XG =XG/L 

MASS =MASS/(0.5*RHO*L**3) 
rz =I1Z/(0.5*RHO*L**5) 


YRDOT =-—0 00176 
YVDOT =-0.03430 
YR =+0.01187 
YV =-0.03896 
YDRS =+0.02345 
YDRB =+0.02345 
NRDOT =-0.00047 
NVDOT =-0.00178 
NR ==0- 01022 
NV =-0.00769 
NDRS #=-0.337*0.02345 
NDRB =+0.283*0.02345 


DH =(1IZ-NRDOT) *(MASS-YVDOT)- 

(MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
AA11=((IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DH 
AA12=((IZ-NRDOT) *(-MASS+YR) - 

(MASS*XG-YRDOT) * (-MASS*XG4NR) ) /DH 
AA21=( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DH 
AA22=((MASS-YVDOT) * (—MASS*XG+4NR) - 

(MASS*XG-NVDOT) *(-MASS+YR) ) /DH 
BB11=((1IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS) /DH 
BB12=((1IZ-NRDOT) *YDRB-(MASS*XG-YRDOT) *NDRB) /DH 
BB21=((MASS-YVDOT) *NDRS-(MASS*XG-NVDOT) *YDRS) /DH 
BB22=((MASS~YVDOT) *NDRB-(MASS*XG-NVDOT) *YDRB) /DH 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


Wiehe *> *) = BNteRe te, fin De 
66 


ILE=1.0/TC 

91=3.0*POLE 

92=3.0*POLE**2 

93=POLE**3 

l=BB1 

l=BB2 

f=—-AD1-(AA11+AA22) 
2=(BB1*AA22-BB2*AA12) 
2=(BB2*AA11-BB1*AA21) 
1=AD3/(BB2*AA11-BB1*AA21 ) 
J=AD2-(AA11*AA22-AA12*AA21)+BB2*K1 
—2=(C1*B2-C2*B1)/(A1*B2-A2*B1 ) 
3=(C2*A1-C1*A2)/(A1*B2-A2*B1) 


IME=100.0 
r=0.1 
=TIME/DT 
Y=TL/DT 


PSIDOT=R 

VDOT =AA11*V + AA12*R + BB1*DEL 
RDOT =AA21*V + AA22*R + BB2*DEL 
XDOT =COS(PSI)-V*SIN(PSI) 

YDOT =SIN(PSI)+V*COS(PSI) 


PSI=PSI+(PSIDOT*DT) 
V=V+(VDOT*DT) 
R=R+(RDOT*DT) 
X=X+ (XDOT*DT) 
Y=Y+(YDOT*DT) 


ICOUNT=ICOUNT+1 

IF (ICOUNT.NE.IY) GO TO 11 
YCON=Y 

ICOUNT=0 


PSIC=—-ATAN(YCON/D) 
DEL=K1*(PSI-PSIC)+K2*V+K3*R 
DEL=D0*TANH(DEL/D0O ) 


WRITE(20,*) T,Y 
WRITE(30,*) T,Y 
CONTINUE 
END 
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PROGRAM THESNEWT.FOR (TWO POINT TIME DELAY ANALYSIS) 


APPENDIX C 


IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 


DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDRS,NDRB,IZ,MASS, 


& NRDOT,NVDOT 


DIMENSION A(4,4),FV1(4),1V1(4),222(4,4) ,WR(4) ,W1I(4) 


DIMENSION DIST(10,10) 
OPEN (11,FILE='TCN8.DAT’ ,STATUS='NEW’ ) 


Vehicle Parameters: 


& 


& 


& 


U= 5) (0) 

RATIO= 0.0 

WEIGHT=435.0 

IZ =45.0 

L =. 3 

RHO =1.94 

G =32.2 

XG =0.0104 

MASS =WEIGHT/G 

YRDOT =-0.00178*0.5*RHO*L**4 

YVDOT =-0.03430*0.5*RHO*L**3 

YR =+0.01187*0.5*RHO*L**3 

YV =-0.03896*0.5*RHO*L**2 

YDRS =+0.02345*0.5*RHO*L**2 

YDRB =+0.02345*0.5*RHO*L**2 

NRDOT =-0.00047*0.5*RHO*L**5 

NVDOT =-0.00178*0.5*RHO*L**4 

NR =-0.01022*0.5*RHO*L**4 

NV =—-0.00769*0.5*RHO*L** 3 

NDRS =-0.337*0.02345*0.5*RHO*L** 3 

NDRB =+0.283*0.02345*0.5*RHO*L**3 

DH =(IZ—-NRDOT) *(MASS-YVDOT)- 
(MASS*XG-YRDOT) * (MASS*XG-NVDOT ) 


AA11=((1IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DH 
AA12=((1Z-NRDOT) *(-MASS+YR)-— 
(MASS*XG-YRDOT) * (-MASS*XG+NR) ) /DH 
AA21=((MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DH 
AA22=((MASS-YVDOT) *(—-MASS*XG+NR)- 

(MASS *XG-NVDOT) *(-MASS+YR) ) /DH 
BB11=((1IZ—-NRDOT) *YDRS—(MASS*XG-YRDOT) *NDRS) /DH 
BB12=((1Z-NRDOT) * YDRB-(MASS*XG-YRDOT) *NDRB) /DH 
BB21=((MASS-YVDOT) *NDRS-(MASS*XG-—NVDOT) *YDRS) /DH 
BB22=((MASS-YVDOT) *NDRB-(MASS*XG-NVDOT) *YDRB) /DH 


BB1=BB11+RATIO*BB12 
BB2=BB21+RATIO*BB22 


WRITE (*,1005) 
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READ Vo, IL 


‘RITE (*,1006) 
EAD (*,*) TC 


PS=1.D-10 
TMAX=100000 


L=0.01 
O 4 I=1,100 
LD = TL * L/U 


DO 1 J=1,100 

TCD=TC*L/U 

POLEC=1.0/TCD 

AD1=3.0*POLEC 

AD2=3.0*POLEC**2 

AD3=POLEC**3 

Al=BB1*U*U 
B1=BB2*U*U 
Cl=-AD1-(AA11+AA22)*U 
A2=(BB1*AA22-BB2*AA12) *U**3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
K1=AD3/( (BB2*AA11-BB1*AA21 ) *U**3) 
C2=AD2-(AA11*AA22-AA12*AA21 ) *U**24+BB2*U*U*K1 
K2=(C1*B2-C2*B1)/(A1*B2-A2*B1 ) 
K3=(C2*A1-C1*A2) /(A1*B2-A2*B1 ) 


Al= (AA11*BB2-AA21*BB1 ) *U*U*U*K1 


A2= (AA11*AA22-AA12*AA21 ) *U*U+(AA11*BB2-AA21*BB1 ) *U*U*U*K3+ 


(AA22*BB1-AA12*BB2 ) *U*U*U*K2-BB2*U*U*K1 
A3=-(AA11+AA22) *U-(BB1*K2+BB2*K3) *U*U 
BO= (AA11*BB2-AA21*BB1 ) *U*U*U*U*K1 
Bl=—-(AA12*BB2-AA22*BB1 ) *U*U*U*K1-BB2*U*U*U*K1 
B2=-BB1*U*U*K1 


AQ=B2*A3-Bl 
BQ=B1*A2-B2*A1-B0*A3 

CQ=B0*A1 

DET=BQ** 2-4 .D0*AQ*CQ 

fF (DET.LT.0.D0) GO TO 4 
ZT1=(-BQ+DSQRT(DET))/(2.0D0*AQ) 
ZT2=(-BQ-DSQRT(DET) )/(2.0D0*AQ) 
IF (ZT1.GE.0.D0) OM=DSQRT(ZT1) 
IF (ZT2.GE.0.D0) OM=DSQRT(ZT2) 
ALPHA1=AQ 

ALPHA2=BQ 

ALPHA3=CQ 

BETA1=-B2 

BETA2=B0+B2*A2-B1*A3 
BETA3=B1*A1-B0*A2 


OMOLD=OM 


BETA=BETA1 *OMOLD**5+BETA2*OMOLD* * 3+BETA3 *OMOLD 
ALPHA=ALPHAI1 *OMOLD* * 4+ALPHA2 *OMOLD* *2+ALPHA3 
BEP=5.D0*BETA1 *OMOLD* * 4+3.D0*BETA2*OMOLD**2+BETA3 
ALP=4.D0*ALPHAI1 *OMOLD* * 3+2.D0*ALPHA2 *OMOLD 
F1=BETA*DSIN(2.D0*OMOLD*TLD)-2.D0*BETA*DSIN(OMOLD*TLD) 

| F2=ALPHA*DCOS(2.D0*OMOLD*TLD)-2.D0*ALPHA*DCOS(OMOLD*TLD) 
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FP1=BEP*DSIN(2.D0*OMOLD*TLD)+BETA*2.D0*TLD*DCOS(2.D0*OMOLD*TLD) 
FP2=-2.D0*BEP*DSIN(OMOLD*TLD)-2.D0*BETA*TLD*DCOS (OMOLD*TLD) 
FP3=ALP*DCOS(2.D0*OMOLD*TLD )-ALPHA*2.D0*TLD*DSIN(2.D0*OMOLD*TLD) 
FP4=-2 .D0*ALP*DCOS (OMOLD*TLD)+2.D0*TLD*ALPHA*DSIN(OMOLD*TLD) 
FPRIME= FP1+FP2+FP3+FP4 

IF (FPRIME.EQ.0.D0) STOP 1112 

IF (F.EO.0.D0) GO TO 2 

OMNEW=OMOLD-F/FPRIME 

OMDIF=DABS (OMNEW-OMOLD) 

IF (OMDIF.LE.EPS) GO TO 2 

OMOLD=OMNEW 

IT=IT+1 

IF (IT.GT.ITMAX) STOP 1111 

GO TO 3 


OM=OMNEW 
XDNUM1=DSQRT( (BO-B2*OM**2) **2+(B1*OM) **2) 
XDNUM2=DSQRT((2.D0*DCOS(OM*TLD)-DCOS(2.D0*OM*TLD) )**2+ 


& (DSIN(2.D0*OM*TLD)-2.D0*DSIN(OM*TLD) )**2) 


XDNUM=XDNUM1 *XDNUM2 
XDDEN=OM* DSQRT( (A1-A3*OM**2) **2+(A2*OM-OM**3) **2) 
XD=XDNUM/XDDEN 
X=XD/L 
DIST(I,J)=X 
WRITE(11,*) X,TL 
TC=TC + 0.015 


1 CONTINUE 

THETL + 9.02 

CONTINUE 

WRITE(11,20) ((DIST(1,J) ,J=1,10) 7 1=1 20) 
FORMAT(10(2X,F4.2) ) 

FORMAT (’ ENTER TIME LAG’ ) 

FORMAT (’ ENTER TIME CONSTANT’ ) 

END 
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